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Abstract 

Wc study the oscillatory dynamics in the generic three-species rock-paper-scissors games with mutations. In 
the mean-field limit, different behaviors are found: (a) for high mutation rate, there is a stable interior fixed 
point with coexistence of all species; (b) for low mutation rates, there is a region of the parameter space 
characterized by a limit cycle resulting from a Hopf bifurcation; (c) in the absence of mutations, there is a region 
where heteroclinic cycles yield oscillations of large amplitude (not robust against noise). After a discussion on 
the main properties of the mean-field dynamics, we investigate the stochastic version of the model within an 
individual-based formulation. Demographic fiuctuations are therefore naturally accounted and their effects 
are studied using a diffusion theory complemented by numerical simulations. It is thus shown that persistent 
erratic oscillations (quasi-cycles) of large amplitude emerge from a noise-induced resonance phenomenon. We 
also analytically and numerically compute the average escape time necessary to reach a (quasi-)cycle on which 
the system oscillates at a given amplitude. 

Keywords: Population Dynamics; Fluctuations and Demographic Noise; Evolutionary Game Theory; Co- 
evolution, Cyclic Dominance and Biodiversity; Cycles and Quasi-Cycles. 
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1 Introduction 



Understanding the mechanisms allow i ng co -e volution and t h e mainte n ance o f biodiversity i s a key i ssue i n 



theoretical biology and ecology (IMavl (|l974l ): 



MichodI ( 2000 ) 



Hubbelll ( 2001 ) 



HakenI (|2002l ): 



Murravl (|2002l ): 



Neall (|2004l )). In this co ntext, cyclic dominance has bee n identified a s a po t ential mech a nism t hat helps pro 



19821) : 



mote species diversity (jPurrett and LevinI (|l994 



1997 ) 



Gilg et al 



(|200lh : 



Kerr et al 



()200% 



( 2002 ) ) and is naturally invest i gated i n the fr a mewo r k of evolutionary 



Hofbauer and Sigmundl ( 1998 ) 



Nowad (|2006l ): 



Szabo and Fath 



game theory (EGT) 



Czaran et al 



Mavnard Smith 



20071 )). The emergence of oscillatory 



(or quasi-oscillatory) behavior is o ne of the most appea l ing and debated phenomena that of t en character i zes co 



evolution in popu l ation d ynamics ( Sineryo and Livelv 



Kirkup and RilevI (|200J); 



Dawkins 



19761) : 



1996): 



BastoUa et al 



Zamudio and Sinervo ( 200C ) 



2OO2I ): 



Hauert et al 



mm; 



Kerr et a 



Mobilia et al 



Oscillatory dynamics has notab l y been ob s erved in predator-prey and host-pathoeei i systems ([Anderson and Mav 



(1986 



19911); 



BerrvmanI (|2002l) 



20021) : 



20071) 'I 



TurchinI (|2003[ )). as well as in genetic networks ([Elowitz and Leibleil (|200ol )). 
Here, we study the oscillatory dynamics of rock-paper-scissors games with mutations and demonstrate that this 
favors the long-lasting co-evolution of all specie s. Rock-paper-scissors games (RPS) -in which r ock crushes scis 



sors, s cissors cut paper, and paper wraps rock (jHofbauer and Sigmundl (|1998I ): 



Nowad (|2006l ): 



Szabo and Fath 



(|2007l )')- ha ve emerg e d as 



ecosystems (Tainakal (119941): 



p aradigmatic mathematica l models in EGT to describ e the cyclic competi t ion in 



Mav and Leonard! (119751): 



Szabo and Szolnoki (2002): 



Reichenbach et al. 



(2006 



2007a . 



2008b ) 



Perc et al 



20071) : 



Claussen and Traulsen 



20081) : 



Peltomaki and Alava 



20081) : 



Berr et al 



mom ) 



In addition to their theoretical relevance, and in spite of their apparent simplicity, it has been argued that the 



the cy clic dominance observed in some communities of lizards (jSinervo and Lively! (|1996I ): 



Zamudio and Sinervo 


li { 


Kerr et al. 




2002 


)). 



In this case, it has been found that cyclic dominance yields coexistence of all species in a spatial setting, 
while in a well-mixed (homogeneous) environment two species go extinct after a short transient. Naturally 
therefore, the mathematical properties of the RPS model and of its variants have recently attracted much 
interest. In an EGT setting, the RPS dynamics is classically described in terms of the rep li cator equation s 



(REs) that is a s et of deterministic rate equations (see below) (jHofbauer and Sigmundl (jl998[ ): iNowakI (|2006l ): 



Szabo and FathI (|2007l) ). The latter essentially predict two types of behavior for the RPS dynamics: there is 
either the stable coexistence of all species, or oscillations of large amplitudes with each species almost taking 
over the entire system in turn and then suddenly almost going extinct. As population dynamics always involves 
a finite (yet, often large) number of discrete entities, it is necess ary to take the effec t s of n o ise into account 



EwensI (|20M)). 



and understand its influence on the ensuing nonlinear dynamics (jNisbet and GurnevI (jl982f ): 
In the context of population models, demographic stochasticity (i.e. intrinsic noise) is naturally accounted by 
adopting an 'individual-based' modeling. The stochastic dynamics is thus implemented in terms of random 
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birth and death events, hke m the Moran process ([MorarJ (|l958l ): iNowakI (j2006[ )). Recently, individual-based 
modeling has been extensively used to study the stochastic RPS dynamics. In fact, a large body of research 
has been concerne d with sp a tially -e xtended systems and the i nfl uence o f the species' movement on t he emerg- 



ing noisy patterns (jTainakal (|l994l ): 



Reichenbach et al 



( 2007al b 



2008b ) 



Peltomaki and Alaval |2008l )). On the 



other hand, for systems with homogeneous (well-mixed) populations, it is well established that intrinsic noise 
drastically alters the repli cator dynamics and recen t studies have focused on the mean time and prob ability of 



extinction of two species ([Reichenbach et al 



Claussen and TraulsenI ( 20081) 



Berr et al 



mom ) 



In addition to selection and reproduction, a third basic evolutionary mechanism is mutation. Th e lat- 



ter is generally regar ded as providing the advantageous traits that survive and multiply in offsprings (jNeal 



(1200J; 



MichodI (l2000r )). Also, from a behavioral perspective, it is recognized that agents never behave perfectly 



ratio nally and it is sens i ble to assume that indi viduals can switch their strategies and therefore undergo muta- 



tions (jAntal et all (|2009 



els considered in Refs. 



Tainaka 



Szabo and FathI (120071) '). While selectioii and reproduction are 



19941) :ISzab6 and Szolnokil (120021): 



Reichenbach et al 



present in th e RPS mod- 



(2006 



2007a) 



Perc et al 



2007); 



Claussen and TraulsenI ( 2008 ) 



Reichenbach et al 



2007b 



2008bl) : 



Peltomaki and Alava 



2008) 



Berr et al. 



20091 )). the latter do not account for mutations. Here, we study the oscillatory dynamics of RPS games with 
mutations and show that the possibility to switch from one strategy (species) to another with a small transition 
rate yields novel oscillatory dynamics and favors long-term coexistence. For population of finite size, we discuss 
the effect of demographic noise on the dynamics and show that it induces persistent (quasi-)cyclic behavior 
characterized by sustained erratic oscillations of non-vanishing amplitude. 

The remainder of the paper is organized along the following lines. In the next section we introduce the generic 
RPS model with mutations (RPSM) and describe its dynamics in the mean-field limit. The (generalized) REs 
are studied in Section 3, where the bifurcation diagram is obtained (Sec. 3.1) and a new region of the parameter 
space characterized by a Hopf bifurcation is identified. The main properties of the resulting limit cycle are 
briefly discussed in Sec. 3.2. Section 4 is dedicated to the stochastic dynamics of the model RPSM in terms of 
an individual-based formulation. In particular, we show that demographic noise can cause quasi-cyclic behavior 
with persistent oscillations of large amplitude (Sec. 4.1), and allow to escape - after a characteristic time that 
is computed - from the coexistence fixed point and reach a given (quasi-)cycle (Sec. 4.2). In the final section. 



we summarize and discuss our results. 



2 Rock-Paper-Scissors with Mutations 

In their essence, all variants of the RPS game aim at describing the co-evolutionary dynamics of three species, 
say A, B and C, in cyclic competition. In this setting, as in the children's game, 'rock' (here species A) crushes 
the 'scissors' (species B), and 'paper' (species C) wraps the 'rock', and 'scissors' cut the 'paper'. We therefore 
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say that A dominates over C, which outcompetes B, which outgrows A and thus closes the cycle. In EGT, 
the interactions are specified in terms o f a payoff matrix V. Generical l y, the cyclic dominance of RPS g ames is 



captured by the following payoff matrix (jHofbauer and Sigmundl (jl998[ ): 
where e > 0: 



Nowakl l 20061) 



Szabo and FathI (|2007t) ). 



V = 



vs 


Rock (A) 


Paper (B) 


Scissors (C) 


Rock (A) 





—e 


1 


Paper (B) 


1 





— e 


Scissors (C) 


— e 


1 






According to this matrix, when a pair of A and B players interacts, the former gets a negative payoff — e while 
the latter gets a payoff 1. In this case, A is dominated by B and its loss is less than B's gain when < e < 1, 
whereas S's gain is higher than A's loss when e > 1. In the same way, C's dominate over _B's and the latter 
prevail against A's, and thus the model exhibits cyclic dominance. In all pairwise interactions, the dominant 
individual gets a payoff 1, whereas the dominated one obtain a negative payoff — e. Therefore, the parameter e 
allows to introduce an asymmetry (when e 7^ 1) in the interactions. When e = 1, one of the player loses what 
the other gains and this perfect balance corresponds to a zero-sum game. 

In EGT, one often considers homogeneous (well- mixed) populations of N individuals, with — > 00. In 
this mean- field limit, the dynamics is usually specified in terms of the REs for the densit i es (or relative 



abund a nces) a(t),b{t) and c(t ) of species A,B and C, respectively (jHofbauer and Sigmundl (|l998t ): 



pood) : 



Nowak 



Szabo and FathI (|2007l )). Introducing the vector s(i), whose components Si, with i € {A,B,C), are 
sa = a(t),SB = b(t) and sc = c{t), the REs read Si = Si[{Vs)i — s.Vs] = Si[7r,; — tt], where the dot stands 
for the time derivative. The important notion of average payoff (per individual) of species i, tt^, has been 
introduced in terms of the payoff matrix as a linear function o f the relative abundances: tt^ = (Vs)i, wherea s 



NowakI pood) : 



7f = S.Vs = - SjiTj denotes the population's mean payoff (IHofbauer and Sigmundl (jl998l) : 
Szabo and FathI (|2007l )). In addition to the above processes of selection/reproduction, we introduce a third 
evolutionary mechanism that allows each individual to mutate from one species to another with rate /x: 



A 



B 



C 



B 



A 



C 



, c 



(1) 



B 



In this setting, the natural generalization of the replicator equations for the model under consideration is 
Si — s^TTi — 7f] + /i(l — 3si), with TTA — c — eb, ttb — a — ec, ttc = b — ea and tt = (1 — e){ab + be + ac). The 
latter are studied in the next sections. 
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Figure 1: (Color online). Typical plots of the time-dependent density (relative abundance) a{t) of species 
A, obtained from Eqs. (|2]) for e = 1.225 and different values of ^ (with a(0) = 0.25,6(0) = 0.40). Here, 
He = 0.0125 (see text). Left: In the absence of mutation rate (/i = 0,A = 0.0375), the fixed point s* is 
unstable and the density (brown, dashed) jumps from to 1. For very low values of the mutation rate, here 
jjL = 0.001 (A = 0.0345), s* is also unstable and a{t) oscillates regularly (magenta, solid) with an amplitude 
that approaches the extreme values and 1, but with a much shorter period than in the case /i = 0. Right: 
For /i — 0.01 (A = 0.0075), the fixed point s* is still unstable and a{t) oscillates regularly about a* ~ 1/3 (red, 
solid), with a period T « 10. When ji = 0.02 > (A = —0.0225), the densities exhibit exponentially damped 
oscillations and converge towards the fixed point value 1/3 (blue, dashed). 



3 The Rate Equations 

Here, we consider the dynamics of the model RPSM in the mean-field limit, where N oo. When the 
population is well-mixed (homogeneous) and every pair of random individuals has the same probability to 
interact, demographic noise can be neglected and the system's dynamics is aptly described by the rate equations 

d = a [c — e6 — (1 — e) {ab + bc + ac}] + — 3a) 

b = b[a- ec- [l - e){ab + bc + ac}]+ ^{l -ib) (2) 
c = c [6 — ea — (1 — e) {ab + 5c + ac}] + /i(l — 3c) 



As the system is comprised of three species, the sum of the population density is conserved, i.e. a(t)+fe(i)+c(i) = 
1. Clearly, this constant of motion allows, say, to set c{t) = 1 — a{t) — b{t) and reduces ([2]) to a system with 
only two variables (here, a{t) and bit)). Thereafter, the properties of these equations are studied and it will be 
shown that mutations yield a new region in the parameter space characterized by a stable limit cycle. 

It is known that in the absence of mutations (/i = 0) the REs admit one interior fixed point, s* = (a*, 6*, c*) = 
(1/3, 1/3, 1/3), and three abs orbing states ((1,0,0), (0,1,0) and (0, 0, 1)). The resulting cyc l ic dyn amics gives 



rise to three kinds of behavior (jHofbauer and Sigmundl (|l998l ): 



NowaklpOOd ): 



Szabo and FathI pOOTl) ): (i) When 



e < 1, s* is the system's only attractor, it is globally stable and the trajectories in the phase portrait spiral 
towards it (s* is a focus), (ii) When e > 1, the interior rest point s* becomes unstable and the flows in the phase 
portrait form a heteroclinic cycle connecting each of the absorbing fixed points (saddles) at the boundary of the 
phase portrait. Clearly, as any fluctuations cause the rapid extinction of two species, the resulting oscillations 
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l<0 

(- Stable interior fixed point 
^ ' (Flows spiralling inwards) 



Hopf bifurcation 




Hc=(e-1)/1^ 



fi^\ Stable oscillations 
^ ' Robust cycles 
O 

1 ? e 

"""(c) |J,=0, 8> 1 : heteroclinic cycles 

Figure 2: (Color online). Bifurcation diagram in the parameter space of the (deterministic) model RPSM. The 
regions (a) and (b) are separated by the critical mutation rate fic = {e — 1)/18 > 0. In region (a), fi > fXc for 
e > 1 and > for e < 1, s* is a stable focus. Below the line fJ-d^), in region (b) where < /i < /ic, the 
system undergoes a Hopf bifurcation (at = /ic) and the dynamics is characterized by a limit cycle and stable 
oscillations. In region (c), where ^ = (no mutations) and e > 1 (dashed line), the REs ([2]) yield heteroclinic 
cycles along with periodic oscillations that are unstable against demographic fluctuations. 



arc non-robust (jMav and LeonardI (|l975l ): iDurrett and LevinI (jl994[ )). (iii) When e ~ 1, tt vanishes (zero-sum 
game) and the interior rest point s* is marginally stable (center). In this case, the quantity a{t)b{t)c(t) is a 
constant of motion and the trajectories in the p hase portrait form c losed orbits around s* , set by the initial 



conditions and not robust against noise (see e.g. (jReichenbach et al 



( 20061) )) 



In the presence of mutations, solving Eqs. ([2]) with a = 6 = c = 0, one finds that s* = (1/3, 1/3, 1/3) is 
the only (interior) fixed point of the system. The absence of absorbing points when /.j > indicates that the 
model under consideration does not yield heteroclinic cycles. In fact, as shown below, in the presence of small 
mutation rate the heteroclinic cycles are replaced by stable cycles resulting from a Hopf bifurcation (region (b) 
in Fig. [2]) and yield persistent oscillatory dynamics of all species, which is a feature of biological relevance. 

3.1 Linear stability analysis 

To investigate the properties of REs (l2|), it is useful to introduce the variables x = {xa, xb)'- 



XA=o,^a*=a and .tr=6 — &*=6 , 

3 3 



(3) 



which measure the deviations from the interior fixed point s* . We notice that xc = c — c* = —{xa + xb)- As 
a first step to gain some insight into the dynamics described by it is useful to perform a linear stability 
analysis. To linear order in terms of a; = [xa, xb)), the REs ^ can be rewritten as ac = A{s*)x, where A{s*) 
is the Jacobian matrix at s* whose (complex conjugate) eigenvalues which are A± = A ± iujQ, with 



A = — 3/i and ujq = 

6 ^ 2x/3 



(4) 
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General results therefore guarantee that s* is (globally) stable when A < 0, i.e. for 18jU > e — 1, whereas s* is 
unstable when A > and a Hopf bifurcation occurs at the critical value fi — fid^) and e > 1 is kept fixed, with 



18 



>0, 



(5) 



The bifurcation diagram of the model is thus summarized in Fig. [51 Such a diagram clearly illustrates the 
effect of mutations and is characterized by three regions. In region (b), where < fi < fic, s* is unstable and 
there is a Hopf bifurcation at /i = /ic (as demonstrated below) which leads to the coexistence of all species 
with stable oscillations of their densities (right panel of Fig. [1] red/solid). On the other hand, in region (a), 

{He ! for e > 1 
, i.e. when A < 0, the densities are exponentially damped, with an amplitude 
, for e < 1 

Aamp(0 ~ e"l^l* that vanishes near s* (right panel of Fig. [1] blue/dashed). In the absence of mutations and 
when e > 1 (region (c) in Fig. [5]), one recovers the scenario characterized by heteroclinic cycles (left panel of 
Fig. [Tl brown/dashed). In this case, the density of each species oscillates, with a large period, between the 
extreme values (absence of a given species) and 1 (presence of only one species). These oscillations are not 
robust against demographic noise: chance fluctuations will unavoidably and quickly drive the system into an 
absorbing state. When A = 0, the fixed point s* is a center and its stability is determined by nonlinear terms. 
Below, we shall consider the dynamics to third order about s* and sho w that it is stable when A = and /i > 0. 



Only in the marginal case A = /i = 0, s* is marginally stable (see e.g. (jReichenbach et al 



mom 



3.2 Hopf bifurcation and limit cycle 

As qualitatively discussed above, at mean-field level, the most interesting effect of mutations is the emergence 
of robust oscillatory behavior (stable limit cycle) resulting from a Hopf bifurcation occurring at low mutation 
rates (right panel of Fig. [TJ red/solid). To study the properties of the resulting limit cycle and take advantage 
of the system's symmetry around s*, it is convenient to perform the linear transformation x — > Sx = y, where 



S = 



^ 

" 2 



As well known from bifurcation theory (jWiggind (|l99Cll ): 



Grimshawl (|l994l) ). one then Taylor 



expands the REs (expressed in the variables y = Sx) to third-order around s* and recasts the resulting rate 
equations into the Hopf bifurcation normal form. Hence, in polar coordinates (r, 9), the normal form of the REs 
(truncated to third order) reads 

f = r(A + ^r^) 



Wo — car 



(6) 
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r^(e,^i) 
0.25^ 




Figure 3: ( Color online). Plot of Vryo = ^A/|/3| > for /i > /ic as function of e and \x according to (|4]) and ([7]). 
We notice that Too increases when e is raised and /i is lowered. 



where 



18a;o(l + 2V3a;o) 



7(1 + e2) + e(13 - 9^) + 9/i(l + 9/^) 



(7) 



/3 = 1 - e - 



6A(1 + 2eV3wo) 



7(1 + £2) + e(13 - 9m) + 9m(1 + 9^) 



One readily recognizes that ([6]) yield a supercritical Hopf bifurcation when A > (for < fi < fic{e)). In fact, 
the parameter /? is the f i rst Lyapunov coefficient a.nd is negative w hen A > and the mutation rate /x is small 



(i.e. n/e < 1) (jWiggind (|l990f ) 



Grimshawl |l994l ): 



strogay diggj)) 



It is insightful to solve the radial component of the Eq. ([5]), whose solution is r{t) 



r(0) 



Thus, 



Vl-r(0)|(e2At_i)-^ 

when A > 0, the interior fixed point s* is unstable and the long-time behavior is r(t) ~ r^o ^1 — — | ^^^^WPP }) ' 
where 



(8) 



The dependence of r^c on the parameters e and fj, is illustrated in Figure [3l where it is shown that Too increases 
when e is raised and /i lowered. Thus, it follows from Eqs. (|6l8p that the REs ^ indeed yield a stable limit 
cycle d'{t) ~ (a(t), b{t)) (with c(t) = 1 — a{t) — b{t)) of radius Too, period T and frequency w, where 



27r Aa 
T= — , with cj = - — 
w L3 



(9) 



We notice that the linear terms contribute to the frequency uj through cjq (natural frequency), while nonlin- 
earity gives rise to an additional contribution (— Aa/|/3|). According to ([6]), the limit cycle o'(t) is approached 



exponentially fast in time, r(t) 
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x,(t) 

Figure 4: (Color online). With Xji{t) = a{t) — ^ and XB{t) = b{t) — i, comparison of the hmit cycle obtained 
from the numerical solution of ^ (red/thick, here 550 < i < 1000) with the predictions of ([5]) [black/thin]. 
The parameters and initial conditions are e = 1.2, ji — 0.01,x'a(0) = 0.025, ^^(O) — 0.040. here A ~ 1/300. 



As they result from a third-order expansion, the predictions of ([6]) are accurate for mutation rates close to 
the critical value /ic (with e fixed) , i.e. for "small" positive values of A and Too as illustrated in Fig. |4l 

To conclude this section, we also consider the case where A < and notice that Eqs. ([5]) thus yield r{t) 
g-|A|t^ Xhis implies that s* is approached in an oscillatory manner with an exponentially damped amplitude, 
namely XA{t) oc e"''^'* (cos (wot) — sin (cL>ot)/\/3) and XB{t) oc e"'^'* sin(a;ot). Furthermore, in the marginal 
case A = (with > 0), s* is still stable but approached in a slow oscillatory manner. In fact, it follows from 
^ that in this case r(t) ^ which implies that the oscillations are damped algebraically as f-^^^. Indeed, 



for A = one finds XA{t) 



cos (ujot) — sin (cjoO /{6\/Jd) and XB{t) sin (ajoi)/(3\/3/it). 



4 The Stochastic NonUnear Dynamics 



We so far have focused on the deterministic description of the model in the mean-field limit. However, as 
virtually all real systems arc of finite size, N < oo, and made of discrete entities, they are influenced by 
stochastic fluctuations. Here, adopting an individual-based approach, we discuss how demographic noise alters 
the deterministic properties of the model RPSM. 

Before mathematically studying the model RPSM in the presence of demographic noise, it is worth gaining 
some insight into its stochastic dynamics. As explained below, the latte r has been im plemented by a birth- 



death process simulated according to the Gillespie algorithm (jGillespid ()1976l Il977l )l. Some characteristic 
results are reported in Figs. [SJ71 The behavior corresponding to the regime (a) of the parameter space (see 
Fig. [2]) is illus trated in Fig. \5\ where dem ographic noise causes erratic oscillations forming "(phase-forgetting) 
quasi-cycles" (jNisbet and GurnevI (|l982l) ) in the phase portrait (sec also Fig. [3 cyan/light grey trajectory). 
When the population size N is large and sample-averaged over many replicates, the amplitude of the stochastic 



oscillations decreases and the predictions of the rate equations arc approached (thick curves in Fig. [5]), but fully 
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Figure 5: (Color online). Erratic oscillations and stochastic dynamics of the model RPSM in the region (a) 
of the parameter space, with e = 0.5, /i = 0.01 and N = 300 (initially Na = 135 and Nb = 66). In single 
realizations, the time-dependent densities (relative abundances) of species A (thin, red/dark grey) and B (thin, 
cyan/light grey) are characterized by persistent stochastic oscillations resulting from a resonance amplification 
occurring at frequency fi* ~ wq ~ 0.43 (see text). When they are sample-averaged (here over 200 replicates), 
the mean densities (solid and dashed thick curves for species A and B, respectively) display damped oscillations 
approaching the value 1/3, in agreement with ([2]). 



recovered only in the mean-field limit N — > oo. The stochastic dynamics in the region (b) of the parameter 
space, where s* is unstable, is illustrated in Fig.[6]and is characterized by persistent oscillations with erratic, but 
non- vanishing, amplitude and phase. As a result, instead of a perfectly closed orbit the trajectories in the pha se 
portrait form a perturbed limit cycle, also called "phase-remembering quasi-cycle" (jNisbet and GurnevI (|l982l )). 
as shown in Fig. [7] (red/dark grey trajectory). Figs. [S] and [H] show that demographic noise affects considerably 
each single realization in the region (a) of the parameter space, where the predictions of the REs are replaced 
by oscillations of non- vanishing amplitude, while it perturbs the cyclic behavior in region (b). Below, we show 
that the "qu asi-cycles" arising in re g ion (a ) of Fig. [5] at low mutation rate result from a resonance amplification 



mechanism (jMcKane and Newman! (|2005h ). We will also see how noise affects the system's attractors that 
can be regarded a s min i ma of an effective p otential well from wh ic h it takes an e norm ous amount of time to 



escape (jGardinerl (|2002l ): 



Kuboetal 



(|l973h 



Dvkman et al. 



(|l994h 



Volovik et al 



iaoogl)) 



To describe the sto chastic nonlinear dynamics of the system, we adopt an individual-based approach and 



consider an urn model (jFelleii (jl968l )) comprising a total of = Na + Nb + Nc individuals, Na are of species 
A, Nb of species B, and Nc of species C. Up to fiuctuation of order 0[N^^/'^) Q (see below), the densities of 
species A, B and C are respectively a = Na/N, b — Nb/N and c ~ Nq/N ~ 1 — a — b. The stochastic dynamics 



is here implement ed according t o the frequency-depen d ent M oran process (fMP) (jNowak et al 



POOGI )) (see also (|Moranl (|l958h : 



|2004h : 



Nowak 



Blvthe and McKand (|2007l )')'). The latter is a Markovian process in which. 



at each time step, two individuals - say one of species A and another of species B - are randomly selected 
and then reproduce according to their fitness, here proportional to the individual's average payoff. Once the 
density-dependent average payoffs are determined, one of the individual, say A, reproduces with a rate T^~^^ 



^In fact, the number of individuals of species A is related to a by Na = o,N + 0{\/N) and similarly with the other species. 
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Figure 6: ( Color online). Stochastic dynamics of the model RPSM in the region (b) of the parameter space with 
e = 1.25, n = 0.008 and N = 600 (initially Na = Nb ^ 240). We notice that the averaged (over 250 rephcates) 
time-dependent densities of species A (thick, solid) and B (thick, dashed) are characterized by oscillations of 
non-vanishing amplitude. The erratic trajectories in red/dark grey and cyan/light grey are single realizations 
(no sample averaging) of the density of species A and B, respectively. 



and the offspring replaces a randomly chosen individual of the other species (B in this example). The fMP 
therefore conserves the total size N of the population and can be regarded as a random-walk with hopping 
rates depending on each species average payoff (fitness). While there are various ways of implementing the 
dynamics of the fMP, we here consider that in the absence of mutations and large (yet finite) population size 
the transition T^^^ from B to A, say, is given by [1 + [tta — tt)] oh. This expression comprises a contribution 
proportional to the average payoff differences {t^a — t^) , t hat accounts f or sele c tion, supplem ented by a constant 



(set to 1) accounting for the background random noise (jNowak et al 



(120041) 



Nowad (|2006l) ). The muhiplying 



factor ab encodes the probability (when N ^ 1) that two individuals of species A and B interact. The effect 
of mutations is thus taken into account by adding a linear term, which finally leads to the following transition 
rate: T^^^{s) = (1 + {tta ~ tt}) ah + [ih. More generally, we consider the following transition rates: 



(1 + {tTj - Tt}) SjSj + ^Sj, 



(10) 



with i 7^ ) G (j4, _B, C) and tt^i = c — th, ttb ~ a — ec, ttc — b — ea and tt = (1 — e){ab + be + ac). 

The stochastic description of the system is encoded in the probability P{Na, Nslt) of having Na and Nb 
individuals of species A and B, respectively , at time t . The q uantity P(Na , Nb', t) obeys the master equation 



associated with the above fMP (see, e.g., (jGardineij (|2002l ): lEwensI (|2004r ))) and specified by the transition 



rates PH)) . The Markovian proce ss defined by the transition rates pUj) has been simulated using the Gillespie 



algorithm (|Gillespid (|l976 . 



19771) ) which was initially introduced to simulate chemical systems through their 



"microscopic reactions" . In the vicinity of s* and for large (yet finite) population size, the master equation of 
the above fMP can be aptly described within a generalized diffusion approximation for the probability density 



P{x,t), where x = (xa,xb) and Xj — — i, j e {A,B) are considered continuous variables (about which 
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Figure 7: (Color online). Phase portrait of the model RPSM in the presence of stochastic fluctuations for 
e = 1.25, = 600, with fi = 0.008 (thick, red/dark grey) and ^ = 0.05 (thin, cyan/light grey), initially 
Na = Nb = 240. The trajectories have been obtained from a sample-average over 250 replicates. When 
II = 0.008, the interior fixed point s* is unstable and each trajectory in the phase portrait forms a "phase- 
remembering quasi-cycle" of erratic non- vanishing amplitude and phase orbiting around a. When /i = 0.05, 
the flow spirals towards the stable interior fixed point s* and erratically wanders in its vicinity forming a 
"phase- forgetting quasi-cycle" (see text). 



we expect fluctuations of order N see below). The temporal development of P{xa, XB',t) is thus described 
by a forward Kolmogorov equation ( FKE) , or Fokker-Planck equation, resulting fro m a size-expansion of the 



master equation (see, e.g., (jGardinerl (|2002f ): 



Performing a (linear) van Kampen expansion 



Ewens (2004): 



Nisbet and GurnevI (119821 ): 



van Kampen 



Traulsen et al 



(2005|))). 



19921 )') about the fixed point s* i n the continuum 



limit (assuming th a t N i s large but finite), the standard procedure (see, e.g., (jGardinen (|2002l ) 



Reichenbach et al 



RiskenI (|l989l ): 



(|2006l) )) leads to the following FKE: 



dtP{x,t) = -9,, [x,A,j{s*)P{x,t)] + -6,,(s*)a,,a,^.P(a;,t), 



(11) 



where we have adopted the summation convention on repeated indices i,j G {A, B). In Eq. (jlip the drift terms 
{A{s*)x)i are obtained from the Jacobian matrix A of t he R Es ([2]) at s* (see Sec. 3.1), while the diffusion 
matrix B is defined by (see, e.g., (IClaussen and TraulsenI (|2008l) )) 



BaaIs*) 
Bbb{s*) 
Bab{s*) 



T 



B^Ai 



(s*) + T'^^"{s*) + T'^-^^(s*) + r^-^<- (s*) 



T 



(s*) +T^^^{s*) +T'^^^{s*) +T^^^{s 



4(1 + 3a^) 

9iV 
4(1 + 3^^) 



9N 



= Ss^(s*) = -{r-4^^(s*) + r^^^(s*)} = - 



2( 1 + 3^) 
9iV ' 



(12) 



As well known, the FKE ([TT|) is equivalent to the following set of linear stochastic (Langevin) differential 
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Figure 8: (Color online). Power spectrum P{n) = {\x{^)\'^) of the model RPSM with e = 0.5, /i = 0.01 and 
N = 300 obtained from numerical simulations and by sample averaging over 200 realizations. We notice an 
isolated peak at characteristic frequency ft* ~ 0.43, in full agreement with the analytical prediction (fTS]). 



equations with Gaussian white noise vector t] — {ija, Vb) (|Gardineij (j2002l ) : iRiskenI (|l989( )): 



XA = Aaa{s*) xa+Aab{s*) XB+rjA 
Xb = Aba{s*) xa+Abb{s*) Xb +T]B, 



(13) 



with (y^i) = and covariance matrix {'rii{'t)Vj{t')) = I3ij{s*)6{t — t'), where i,j £ {A, B) and (...) denotes 
the ensemble average over a large number of replicates. As anticipated, it follows from (jl2ll3p that the noise 
strength is oc N^^^^ and yields fluctuations of order 0{N^^^^) around the deterministic values of xa and xb- 
Below, the discussion of stochastic effects on the evolutionary dynamics of the model are centred on two 
aspects: (i) we show how demographic noise leads to the emergence of quasi-cycles; (ii) and we compute the 
average time for a system with the same initial density of each species to "escape" from the interior rest point 
and reach a given cycle. 



4.1 Quasi-cycles and noise-induced resonance amplification 



In this section, we are especially interested in the stochastic dynamics in region (a) of the phase parameter, where 
A < and s* is a stable fixed point. Here, we investigate one of the most intriguing effects of intrinsic noise, 
which yields persistent erratic oscillations around s* (sec Fig.[5|) and therefore considerably alters the predictions 
of the REs In fact, while one could naively expect only small corrections (of order iV~^/^, see Eqs. (|12ll3p ') 
to the deterministic predictions, it has been suggested that demographic noise can be sufficient to perturb 
the stationary state predi cted by the mean-field analysis and to produce persistent erratic oscillations (see 



e.g. (jBartlettl (|l957 



19601 ))). This is indeed illustrated by the numerical simul ations reported in F i g. [51 It 



has recently been shown that such a behavior, often referred to as quasi-cyclic (jNisbet and Gurnev 



ilMl), 
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results from a resonant amplification of demographic (intrinsic) fluctuations (jMcKane and NewmanI (|2005D I. 
To analytically demonstrate this phenomenon, we shall compute the power- spectrum P(f2) of the system. The 
power-spectrum is one of the most useful tools t o look for oscillations in noisy data a nd is related to the 



auto-correlation functions of (stationary) systems (jNisbet and GurnevI (|l982l ): 



Gardineij (|2002l )). Introducing 



x{Vt) = 1^00*^^ ^ *^^*a;(t), the Fourier transform of x{t), the power spectrum is given by the following 
ensemble average: P{n) = (|£(0)p) = 2{\xa{^)\'^) = 2{\xb{^)\'^) ■ Taking the Fourier transform of Eqs. (HH) 
and then averaging the square modulus of x{n) one finds: 



Pin) 

with 9ilQ 



8(1 + 3a^) 
9A^ 



(02 - niY + {2xnY 



1 + 2\/3ewo + 9^(9/^ + 1 - e). 



(14) 



For low mutation rate fi (i.e. when n/e ^ 1), 51q > 2X^ and the power spectrum P{il) has a single peak (for 
n > 0) at the characteristic frequency 



n* 




(15) 



as illustrated in Fig. [51 When the mutation rate fi is low, the value of ft* is very close to (yet different from) luq 
and Qq (gl [H]) y. When (A/fio)^ <C 1, the above expression simplifies to give fi* = fio(l + ©(A^/ilg)). It is also 
worth noticing that at the onset of the Hopf bifurcation, i.e. when A = (with /i > 0), the denominator of P( fl) 
vanishes for O = JIq f^d in this case the peak in th e powe r spect rum is replaced by a pole at = SI* = Qq. 

While on general grounds (law of large numbers ( Fellerl ( 19681 ))) one expects oscillations of the order N'^^"^, 
the presence of a peak in the power spectrum corresponds to a noise-induced amplification due to a resonance 
effect. In fact, there is a resonant behavior when there exists a particular frequency for which the denominator 
of the power spectrum P{n) is small. Here, the denominator of the power spectrum takes its minimal value at 



(16) 



which simplifies when (A/SIq)^ ^ 1 and also gives fires = fio(l + ^(A^/rjQ)). In such a regime, where ft* ~ 
ilj-cs — ^0 , there is a very suggest ive explanation of the noise-i nduced resonant amplification by comparison 



with a simple mechanical system (jMcKane and NewmanI (120051 ) ). Indeed, one can consider a linear damped 



harmonic oscillator of natural frequency flo, with damping constant C, driven by a force oscillating with frequency 
n. One can thus show that such a mechanical system oscillates at frequency SI, with amplitude Aamp(fi) oc 



[(S72 _ j^q) + (C^i)^] and yields a resonant frequency fires = y ^0 ~ Thus, the quantity 2A can be 
interpreted as the damping constant C of the mechanical system. Therefore, provided that 2A is much smaller 

^E.g., for £ = 0.5 and ^ = 0.01, one finds: ujq ~ 0.4330, whereas Q* ~ 0.4328 and Co ~ 0.4476. 
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than Hq (i.e. A^/fig <C 1), the model's dynamics is essentiahy the same as for a hnear damped harmonic 
oscillator with a resonant frequency fl^cs — ^o- Yet, while the driving frequency has to be tuned to achieve 
resonance in the mechanical system, this is not the case for the noisy system that we are considering. In fact, 
as the model RPSM is driven by internal (Gaussian) white noise (demographic stochasticity) which covers all 
frequencies, the resonant frequency of the system will be excited without need of any external tuning. 

Other fruitful quantities to measure repetitiveness and (quasi-)periodicity in a fluctuating population are 
the autocorrelation functions {xAiT + t)xAiT)) = {xsit + T)xB(t)) and (x^(t + t)xB(r)) ~ — (^^(t + t)a;^(T)). 
The latter are related to th e the power spec trum by the Wiener-Khinchin theorem and can be computed by the 
Fourier transform of P{^) (jGardineii (|2002[ )). which yields (for t > r ^ 1) 



{xAir + t)xA{T)) 

{xa{t + t)xB{T)) 



4(1 + 3m) e- 



-|A|t 



3A^ |A| 
2(1 + 3m) e-l^l* 



3A^ 



lAI 



cos (woOi 

sin (woi) - cos (iii;oi)| 



(17) 



Following Ref. (jNisbet and GurnevI (Il982l )). fluctuations arc non-cyclic if the auto-correlation functions de- 
cay monotonically to zero, whereas they lead to quasi-cycles if the auto-correlations oscillate. In the latter 
case, quasi-cycles are "phase-forgetting" if the oscillations of the auto-correlation functions are damped, as in 
Eqs. (|17p . and "phase-remembering" when their amplitude do not vanish. Following this classification, it is clear 
from P7)) and Figs. [5] and [5] that the quasi-cycles in the region (a) of the parameter space are "phase- forgetting" , 
while those in the region (b) are "phase-remembering" . 

An analysis similar to that of this section can be carried out in the region (b) of the phase parameter by 
performing a van Kampen expansion around the limit cyc le a (i). This l eads t o deal with a stochastic version 



of non-autonomous linear differential equations (see, e.g., (jBoland et al 



mm)) 



4.2 Average escape time 

Another intriguing effect of demographic noise is related to the pivotal concepts of attractor and stability. 
As discussed in the previous sections there are two different types of stable attractors in the model RPSM: 
s* = (1/3,1/3,1/3) is the only attractor of the system when A < 0, whereas all trajectories in the phase 
portrait approach the limit cycle & = {a{t),b{t)) when A > 0. This scenario predicting drastically different fate 
depending on whether A < or A > has to be revised when fluctuations are accounted. In this case, within a 
probabilistic setting, attractors can be regarded as minima of a potential well fr om which i t is a l ways possible t o 



escape after some c haracteristic time (wh i ch can be enormously long, see, e.g., (jGardineij (|2002l ) 



Kubo et al 



( 19731) : 



Dvkman et al. 



(Il994h : 



Volovik et al 



RiskenI (|l989l) : 



(|2009l) )). Stochasticity affects the concept of stability 



even if the system is initially at the fixed point s* . In fact, due to demographic fluctuations, the species densities 
deviate from s* and yield quasi-cycles resulting in persistent erratic oscillations of non-vanishing amplitude. It 
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Figure 9: (Color online). Average escape time Tesc(-R) to reach the cycles C{R) starting from s* , for R = 
0.10 «»,0.12 (□),0.15 (•),0.17 (x) and large population size N. Left: Plot of In (A^Tcsc(i?)) versus N in the 
case A < 0, with e = 0.5 and /i = 0.01, and sample-average over 100 replicates. We find a linear dependence 
that yields Tosc(^) oc e'-^^'-'^^ ^ /N. Numerical results (filled symbols) are also compared to the predictions of 
Eqs. dUI) (sohd) and ^ (dashed), yielding Ci{R) cx i?^ Right: Plot of exp(2Arosc(-R)) versus N in the case 
A > 0, with e — 1.5 and fi ~ 0.01 and sample-average over 400 replicates. The observed linear dependence 
yields Tf,sc{R) oc ln(C2(-R) A^)/(2A). Numerical results (filled symbols) are also compared to the predictions of 
Eqs. (HD) (solid) and ([^ (dashed lines, almost indistinguishable from the sohd lines), yielding C2(i?) oc R^. 
Inset: Marginal case A = with parameters e = 1.18 and n = 0.01, and sample-average over 100 replicates. 
Plot of R-'^ In (Tesc(i?)) versus N for values of i? = 0.30 (★), 0.35 (A), 0.40 (V), 0.45 {+). For large population 
size, the (approximate) linear dependence In (Tosc(^)) oc 0.23A^ (dashed line as a guide for the eyes) yields 



is therefore biologically relevant to assess the robustness of the population composition and understand how 
noise affects the co-evolution of sub-populations that initially coexist with the same density. Here, to further 
assess the influence of intrinsic fluctuations when N is large yet flnite, we are interested in the average time to 
escape from s* and reach a specifled separating distance from it. In other words, we compute the average time 
Tcsc{R) for a trajectory starting at s* to reach a cycle C{R) on which oscillations around s* are of amplitude 
2R/^/3 (sec below). Below, we discuss our results, that are summarized in Fig. [SI and also present an analytical 
approach - based on the diffusion theory and van Kampen expansion - to compute TcsciR), that is an accurate 
approximation of TcsciR) around s* . 

Using the Gillespie algorithm, we have computed numerically the average escape time Tcsc{R) necessary 
to reach a cycle C{R) starting from the interior fixed point a* = 6* = c* = 1/3 (same initial density of each 
species). To exploit the symmetry of the model (see Sec. 3.2), we work with the variables y ~ {yA,yB) = Sx 
and consider cycles C{R) of parametric equation i/a — Rcoscf), i/b ~ Rs'mcf). In the original a;— variables, the 
equation of C{R) is xa = i? (cos — sin (/)/-\/3) , xb = (2i?/-\/3) sin0 and the amplitude of the oscillations is 
2R/^/3. In our numerical computations, we have considered systems of population size N, with N = 300 — 4000, 
and values of R ranging from 0.10 to 0.45. Each simulation has been sample-averaged over 100 to 400 replicates. 
As we are mainly interested in the linear regime around s* , where analytical calculations are amenable (see 
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below), particular attention has been dedicated to values of i? « 0.10 — 0.20 (i.e. amplitude « 0.11 — 0.23). 
Typical results are reported in Fig. |9l where one essentially distinguishes two cases. If A < 0, s* is stable and 
In (iVTosc(^)) grows linearly with N (for large N) with a slope depending on R, as illustrated in Fig. [9] (left), 
which yields Tcsc{R) oc e'^^'^^^/A^. On the other hand, when s* is unstable (A > 0) and N is large, we find 
that exp(2Arosc(^)) varies linearly with N (slope depending on R), as shown in Fig. [3] (right). This leads to 
Tcsc{R) oc ln(C2(i?) A^)/(2A). In the above asymptotic expressions, Ci{R) and C2{R) are monotonic functions 
of R and comparison with analytical calculations around s* suggests that Ci ^ C2 oc R^ (see below). For the 
marginal case A = and large N , numerical results are reported in the inset of Fig. [S] (left) and we find that 
(InT'osc(^))/^^ (approximately) exhibits a linear dependence on iV (with constant slope). In the case A = 0, 
one thus finds T^sdR) ^ e^^^- ^ , where C3 is a constant. 

In the realm of the diffusion theory and van Kampen expansion, the calculation of the (approximate) average 
escape time %s c{x) can be rega r ded as a first-pa ssage problem to an absorbing cycle C{R) starting from a system 



initially at s* 



Gardiner! (|2002[ ): iRedneij (|200l[ )). To obtain an equation for T^sc it is useful to consid er the so 



called backward Kolmogorov equatio n (BKE), which is the adjoint of Eq. ([TT|) and reads (see, e.g., (| Gardiner 



(I2OO2I) 



RiskenI (|l989f ) 



Redneij (12001^ 1: 



-dtPb{x,t) = CbPb{x,t), 

with £b{x) = Aij{s*)xj d^, + ^B.i.j{s*)dx,dx^. 



(18) 



As in Section 3.2, to reveal the polar properties of the system in the vicinity of s* , it is natural to perform the 

linear transformation a; — > y = Sx and work in the y— variables. Under this transformation, one has Pb{x, t) — > 

/ 



Pb{y,t), A{s*) SA{s*)S-' 



A -ujo 



and B{s*) SB{s*)S'^ 



A 



1+3m 

3N 




( Reichenbach et al. 



(|2006l )). Thus, in the y— variables, the differential operator of Eq. p8| becomes Cb{x) Cb{y) — \{yAdy 



6N y^VA 



dy^). To exploit the system's symmetry around s* , we then adopt 



the polar coordinates {pT<j}), with yA ~ pcoscj) and ys — psiiKp. The radial BKE thus reads dtPb{p,t) 
Cb{p)Pb{p,t), with 



Cb{p) = 



Xp 



l + 'ip\ 1 
P 



6N 



1 + 3/i 

6A^ 



di 



(19) 



To derive the above radial operator, we have assumed that the initial radial symmetry of the probability 
distribution (initially the system is at s*, without any angular dependence) is approximately preserved by the 
dynamics. The above approach is certainly valid at linear order around s* and is consistent with the van 
Kampen linear expansion that has led to PT|) and With the equation for the average escape time 
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Tcscip) at a distance p from s* is given by (jGardineii (|2002D : IRiskenl (|1989f )) 



1 = Chip) Tcscip)- 



(20) 



In the framework of diffusion and van Kampcn approximations, we now analytically compute the mean escape 
time 7^sc(^) for p to attain the value R starting from s* (i.e. initially p = 0), which corresponds to the mean 
time necessary to reach the cycle C{R). Our analytical treatment is valid in the linear regime around s* , where 
^ ~ Wo, and the amplitude of the oscillations on C{R) is 2i?/\/3. In our theoretical setting, the average escape 
time TcsciR) is obtained by sol ving; Eg. (I20p su bject to a reflecting and an absorbing boundary conditions at 



p — and p = R, respectively ([Gardiner 



20021 ) V When A 7^ 0, the solution to this problem can be expressed 



in terms of the function 4'(z) = zexp (^j^^ 



%sc{R) 



6N 
1 + 3/i 



i-j-'Au nil 

— (1-e") 
u 



(21) 



The solid curves in Fig. [5] show the comparison between the prediction (PT|) and numerical results. We notice 
that for small values of R, e.g., R = 0.10 — 0.15, there is an excellent agreement between (pij) and the results 
of numerical simulations (both in the cases A < and A > 0). For larger values of R (e.g. for R > 0.17), 
nonlincarities and non-vanishing angular dependence cause some systematic deviations from the results of 
numerical computations (discrepancies of « 5% — 10% in the results of Fig. [5] for R = 0.17), but (PTjl still 
provides the correct (qualitative) behavior and a reasonable approximation of TcsciR)- the asymptotic limit 
of large population size N, two distinct behaviors can be obtained from Eq. ([?!]) : 

• When A < 0, s* is stable and the main contributio n to ([2T]) is given by T^sc (R) — o pry Ei f ^^rj^/T") ' '^h^re 



Ei(a;) = dt (e* /t) is the exponential integral ( Abramowitz and StegunI ( 1965 )). From the properties 
of this function, we infer the following asymptotic behavior (for |A|iVi?^ ^1): 



TcsciR) 



f l + 3p 
\6iXR)^N 



exp 



3\X\R^N 
1 + 3/x 



(22) 



This result predicts that the mean escape time increases dramatically with the total number of individuals 
N and, when |A|iVi?'^ ^ 1, is proportional to the exponential of the population size. The dashed lines in 
Fig. [S] (left) illustrate that (|22p is a very good approximation of Tcsc iR) in the linear regime around s* , e.g. 
for R = 0.10 — 0.15, when ^ llfj. In particular, we notice that ([22]) coincides with the aforementioned 
asymptotic expression of TcsciR) with Ci = 3|A|i?^/(l + 3/i). For R = 0.17, when N is large, nonlinear 
effects are important and in Fig. |9]thc exponent of (f22|) overestimates that of TcsciR) by ~ 5%. 



3ln Fig. |9] (left), 3\\\R'^N/{1 + 3^) lies between ^ 1.0 (for = 0.1 and N = 300) and Si 9.5 (for _R = 0.17 and = 1000). 
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• When A > 0, demographic fluctuations always cause the departure from s* (unstable), and all trajectories 
in the ph ase portrait are attracted by the limit cycle &. In this case, from the properties of the exponential 



integral (jAbramowitz and Stegun 



|l965l) l. for XNR"^ > 1, the main contribution to ((HI) is 



%sc{R) ^ ^ 



(23) 



where 7em = 0.57721... is Eulcr-Mascheroni's constant. This result predicts that, when XNR^ ^ 1, the 
mean escape time grows logarithmically with the population size N, i.e. Tcsc{R) (ln-/V)/2A, and the 
value of R contributes to subdominant corrections. Here also, the dashed lines in Fig. ^ (right), almost 
indistinguishable from solid curves, agree very well with the numerically computed Tcsc{R) around s* (i.e. 
for R = 0.10 — 0.15). We notice that, when C2 = 3Ai?^/(l + 3/i), the prediction ((231) coincides with the 
asymptotic expression of TcsdR) inferred from Fig. IH] (right). 

It is also worth noticing that results (PT|) - (P5|) are also valid in the absence of mutations (i.e. when fi = Q). 

Besides the above cases, a special situation arises when A = 0, /i > and s* is a center. In fact, in 
Sec. 3.2, we have seen that in this case s* is an attractor, but the dynamics in its vicinity is very slow. To 
account for such a slow dynamics it is not sufficient to consider (jlll llSp - obtained by linearization about 
s* - but nonlinearities ha ve to be accounted. Generalizing the above approach and keeping the cubic terms 



in the drift contribution (jCremeij (|2008l )). leads to (|20|) with the following differential operator [instead of 
p9|]: Chip) = {A + /3p^}p + ^ dp + (^H^) <9p. When A = 0, proceeding as above, one still finds an 

exponential asymptotic behavior (fiR'^N ^ 1): %sc{R) ^ iV~^/^ cxp(27/ii?*Af/(l + 3/i)). This estimate, that is 
consistent with the numerical results reported in Fig. 15] (left, inset), implies that (at given R and N) the mean 
escape time T^sc around s* is significantly shorter in the marginal case than in the case A < 0. 



5 Summary and Conclusion 

This work has been motivated by the importance of further understanding the mechanisms at the origin of per- 
sistent and erratic oscillatory dynamics in population biology, which is frequently observed but whose theoretical 
origin is often debated. Furthermore, this work concerns the effects of mutations in rock-paper-scissors (RPS) 
games, which are regarded as paradigmatic models to describe the co-cvolutionary dynamics of systems with 
codominant interactions. It has recently been proposed that codominance can help promote the maintenance of 
biodiversity, and it is therefore of biological relevance to study generalizations of the RPS model to understand 
which ingredients favor the long-term coexistence of the sub-populations. 

Here, we have investigated the oscillatory dynamics of the generic rock-paper-scissors games with mutations 
in a well-mixed (homogeneous) population of N individuals. Our study has been carried out in the mean-field 
limit {N — > 00) and in the presence of demographic noise {N large but finite). In addition to the regions 
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of the parameter space associated with a stable interior fixed point s* and heterochnic cycles, the possibility 
for the individuals to switch from one strategy (species) to another with a small transition rate yields a new 
oscillatory behavior associated with a limit cycle and resulting from a Hopf bifurcation. In the mean-field limit 
(N — > oo), we have recast the system's rate equations in a normal form and studied the main properties of the 
ensuing limit cycle. When the population size is finite, demographic noise has to be taken into account and the 
dynamics becomes genuinely stochastic and nonlinear. The latter situation has been described in terms of an 
individual-based formulation, with the stochastic dynamics implemented according to a Moran process. The 
influence of intrinsic noise on the dynamics has been analytically assessed within a diffusion theory approxima- 
tion complemented by numerical simulations (using the Gillespie algorithm). We have shown that demographic 
noise transforms damped oscillations into ( "phase- forgetting" ) quasi-cyclcs and perturbs the amplitude and the 
phase of l i mit cycles ( th at become "phase-rememb e ring" quasi-cycles) . As a lso ob s erved in other s ystern s (e.g.. 



(Bartlett (1957 



19601); 



Nisbet and GurnevI (|l982l ): 



McKane and NewmanI (|2005l ) 



Boland et al 



(I2OO8I))), the 



effect of stochasticity is particularly striking in the region of the parameter space where s* is stable. There, 
we have found sustained oscillations driven by demographic stochasticity. In fact, while fluctuations are of 
order N^^^^ , for small mutation rate, the amplitude of the oscillations is amplifled by a resonance amplification 
caused by internal (white) noise. We have shown that this phenomenon translates into an isolated and well- 
marked peak in the power-spectrum and have also computed the autocorrelation functions of the system. To 
further assess the robustness of the co-evolutionary dynamics against noise, we have computed numerically and 
analytically - using a diffusion approximation and a mapping onto a first-passage problem - the average escape 
time necessary to reach a cycle on which the oscillations attain a given amplitude. We have therefore shown 
that such a mean exit time grows either exponentially or logarithmically with the system size, depending on 
whether the interior fixed point is dcterministically stable or unstable, respectively. 

In summary, in the presence of mutations the RPS model yields sustained oscillations in every region 
of the parameter space. The latter are caused by demographic noise and/or result from a Hopf bifurcation. 
Therefore, mutations in the RPS model ensure that all species co-evolve and oscillate in time without ever going 
extinct. One can therefore regard mutations as a mechanism favoring the oscillatory dynamics in communities 
characterized by cyclic co-dominance. 

In future research, it would be interesting to consider spatially-extended versions of this model and study 
how the population will self-organize and which kinds of spatial-temporal patterns would emerge. In particular, 
as it has recently been shown that mobility can aflect co- ev olution in RPS models by mediating between mean- 



fleld and stochastic dynamics (IReichenbach et al 



(2007aUb 



3 



20080 )). it would be relevant to investigate variants 



of the model RPSM where individuals are allowed to move and interact in space. 
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